library(magrittr)
library(tidyverse)
library(Seurat)
library(readxl)
library(cowplot)
library(colorblindr)
library(viridis)
library(magick, lib.loc = "/home/uhlitzf/miniconda3/lib/R/library")
library(ggpubr)
library(nichenetr)
## load global vars: 
source("_src/global_vars.R")

# meta_tbl
# clrs
# markers_v7
# markers_v7_super
# cell_type_super_lookup

names(clrs$cell_type) <- str_replace_all(names(clrs$cell_type), "\\.", " ")
names(clrs$cell_type) <- str_replace_all(names(clrs$cell_type), "Ovarian", "Ov")
seu_cc <- read_rds("/work/shah/uhlitzf/data/SPECTRUM/freeze/v7/Ovarian.cancer.super_processed_filtered_annotated_sample.rds")
seu_tc <- read_rds("/work/shah/uhlitzf/data/SPECTRUM/freeze/v7/T.super_processed_filtered_annotated_sample.rds")
seu_ml <- read_rds("/work/shah/uhlitzf/data/SPECTRUM/freeze/v7/Myeloid.super_processed_filtered_annotated_sample.rds")
seu_fb <- read_rds("/work/shah/uhlitzf/data/SPECTRUM/freeze/v7/Fibroblast.super_processed_filtered_annotated_sample.rds")

seu_ml$cell_type <- "Myeloid.cell"

Idents(seu_cc) <- seu_cc$cluster_label
Idents(seu_tc) <- seu_tc$cluster_label
Idents(seu_ml) <- seu_ml$cluster_label
Idents(seu_fb) <- seu_fb$cluster_label

seu_cc$consensus_signature <- deframe(select(meta_tbl, sample, consensus_signature))[seu_cc$sample]
seu_tc$consensus_signature <- deframe(select(meta_tbl, sample, consensus_signature))[seu_tc$sample]
seu_ml$consensus_signature <- deframe(select(meta_tbl, sample, consensus_signature))[seu_ml$sample]
seu_fb$consensus_signature <- deframe(select(meta_tbl, sample, consensus_signature))[seu_fb$sample]

seu_merge <- seu_cc %>% 
  merge(seu_tc) %>% 
  merge(seu_ml) %>% 
  merge(seu_fb)
ligand_target_matrix <- read_rds(url("https://zenodo.org/record/3260758/files/ligand_target_matrix.rds"))
lr_network <- read_rds(url("https://zenodo.org/record/3260758/files/lr_network.rds"))
weighted_networks <- read_rds(url("https://zenodo.org/record/3260758/files/weighted_networks.rds"))
## one receiver against all senders --------------------------------------

nichenet_output_cell_type <- read_rds("/work/shah/uhlitzf/data/SPECTRUM/freeze/v7/nichenet/nichenet_output_cell_type.rds")
nichenet_output_cluster_label <- read_rds("/work/shah/uhlitzf/data/SPECTRUM/freeze/v7/nichenet/nichenet_output_cluster_label.rds")

## pairwise receiver sender comparisons --------------------------------------

nichenet_output_cell_type_pw <- read_rds("/work/shah/uhlitzf/data/SPECTRUM/freeze/v7/nichenet/nichenet_output_cell_type_pairwise.rds")
# nichenet_output_cluster_label_pw <- read_rds("/work/shah/uhlitzf/data/SPECTRUM/freeze/v7/nichenet/nichenet_output_cluster_label_pairwise.rds")

1 Ligand receptor pairs per condition

extract_df <- function(df_name) {
  
  df1 <- map_depth(nichenet_output_cell_type, 2, df_name) %>% 
    map_depth(1, bind_rows, .id = "comp") %>% 
    bind_rows(.id = "receiver") %>% 
    mutate(interaction_lvl = "cell_type",
           sender = "all")
  
  df2 <- map_depth(nichenet_output_cluster_label, 2, df_name) %>% 
    map_depth(1, bind_rows, .id = "comp") %>% 
    bind_rows(.id = "receiver") %>% 
    mutate(interaction_lvl = "cluster_label",
           sender = "all")
  
  df3 <- map_depth(nichenet_output_cell_type_pw, 3, df_name) %>% 
    map_depth(2, bind_rows, .id = "comp") %>% 
    map_depth(1, bind_rows, .id = "receiver") %>% 
    bind_rows(.id = "sender") %>% 
    mutate(interaction_lvl = "cluster_label")
  
  df_all <- bind_rows(df1, df2, df3) 
  
  return(df_all)
  
}

ligand_activties_tbl <- extract_df("ligand_activities") 
ligand_receptor_tbl <- extract_df("ligand_receptor_df_bonafide") 

ligand_receptor_tbl %>%
  filter(ligand %in% c("CD274"))
receiver comp ligand receptor weight interaction_lvl sender
CD8.T.CXCL13 HRD_Dup_FBI CD274 PDCD1 0.7308595 cluster_label all
CD8.T.ISG HRD_Del_FBI CD274 PDCD1 0.7308595 cluster_label all 
ligand_receptor_tbl %>% 
  filter(weight > 0.5) %>% 
  ggplot(aes(interaction(ligand, receptor), receiver, fill = weight)) + 
  geom_tile() + 
  facet_grid(~comp) +
  scale_fill_viridis()

ligand_receptor_tbl %>% 
  filter(weight > 0.5) %>% 
  group_by(ligand, receptor, receiver) %>% 
  mutate(n = length(receiver)) %>% 
  filter(n == 1) %>% 
  ggplot(aes(receiver, interaction(ligand, receptor), fill = weight)) + 
  geom_tile() + 
  facet_grid(comp~., scales = "free_x", space = "free_x") +
  scale_fill_viridis()

2 Standard nichnet plots

nichnet_plot_wrapper <- function(x) {
  
  p1 <- plot_grid(x$HRD_Dup_FBI$ligand_receptor_heatmap_bonafide + labs(title = "HRD_Dup_FBI"),
                  x$HRD_Dup_HRD_Del$ligand_receptor_heatmap_bonafide + labs(title = "HRD_Dup_HRD_Del"),
                  x$HRD_Del_FBI$ligand_receptor_heatmap_bonafide + labs(title = "HRD_Del_FBI"), ncol = 3)
  
  p2 <- plot_grid(x$HRD_Dup_FBI$ligand_activity_target_heatmap + labs(title = "HRD_Dup_FBI"),
                  x$HRD_Dup_HRD_Del$ligand_activity_target_heatmap + labs(title = "HRD_Dup_HRD_Del"),
                  x$HRD_Del_FBI$ligand_activity_target_heatmap + labs(title = "HRD_Del_FBI"), ncol = 3)
  
  return(plot_grid(p1, p2, ncol = 1))
  
}

2.1 T.cell => Ovarian.cancer.cell

nichnet_plot_wrapper(nichenet_output_cell_type_pw$T.cell$Ovarian.cancer.cell)

2.2 Ovarian.cancer.cell => T.cell

nichnet_plot_wrapper(nichenet_output_cell_type_pw$Ovarian.cancer.cell$T.cell)

2.3 T.cell => Myeloid.cell

#nichnet_plot_wrapper(nichenet_output_cell_type_pw$T.cell$Myeloid.cell)

2.4 Myeloid.cell => T.cell

#nichnet_plot_wrapper(nichenet_output_cell_type_pw$Myeloid.cell$T.cell)

2.5 Myeloid.cell => Ovarian.cancer.cell

#nichnet_plot_wrapper(nichenet_output_cell_type_pw$Myeloid.cell$Ovarian.cancer.cell)

2.6 Ovarian.cancer.cell => Myeloid.cell

#nichnet_plot_wrapper(nichenet_output_cell_type_pw$Ovarian.cancer.cell$Myeloid.cell)

2.7 any => CD8.T.CXCL13

nichenet_output_cluster_label$CD8.T.CXCL13$HRD_Dup_FBI$ligand_activity_target_heatmap + labs(title = "HRD_Dup_FBI")

nichenet_output_cluster_label$CD8.T.CXCL13$HRD_Dup_HRD_Del$ligand_activity_target_heatmap + labs(title = "HRD_Dup_HRD_Del")

3 Interesting observations

get_data_wrapper <- function(seu_obj) as_tibble(FetchData(seu_obj, c("umapharmony_1", "umapharmony_2", "cell_type", "cluster_label", "IL27", "IL6R", "CD274", "PDCD1")))

plot_data <- lapply(list(seu_cc, seu_fb, seu_ml, seu_tc), get_data_wrapper) %>% 
  bind_rows

plot_wrapper <- function(gene) {
  gene <- enquo(gene)
  ggplot(plot_data, aes(umapharmony_1, umapharmony_2, color = !!gene)) +
    geom_point(size = 0.1) + 
    facet_wrap(~cell_type) +
    scale_color_viridis() + 
    theme_void()
}

plot_wrapper(IL27)

plot_wrapper(IL6R)

plot_wrapper(CD274)

plot_wrapper(PDCD1)